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Abstract 

In a semiconductor superlattice with long scattering times, damping of Bloch oscillations due to 
scattering is so small that nonlinearities may compensate it and Bloch oscillations persist even in 
the hydrodynamic regime. To demonstrate this, a Boltzmann-Poisson transport model of miniband 
superlattices with inelastic collisions is proposed and hydrodynamic equations for electron density, 
electric field and the complex amplitude of the Bloch oscillations are derived by singular pertur- 
bation methods. For appropriate parameter ranges, numerical solutions of these equations show 
stable Bloch oscillations with spatially inhomogeneous field, charge, current density and energy 
density profiles. These Bloch oscillations disappear as scattering times become sufficiently short. 
For sufficiently low lattice temperatures, Bloch and Gunn type oscillations mediated by electric 
field, current and energy domains coexist for a range of voltages. For larger lattice temperatures 
(300 K), there are only Bloch oscillations with stationary amplitude and electric field profiles. 

PACS numbers: 72.20Ht, 73.63.-b, 05.45.-a 
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I. INTRODUCTION 



Bloch oscillations (BOs) are coherent oscillations of the position of electrons inside energy 
bands of a crystal under an applied constant electric field —F. Their frequency is propor- 
tional to the field F and to the lattice constant I: weioch = eFl/h. BOs were predicted 
by Zener in 1934 as an immediate consequence of the Bloch theorem [l], but they were 
not experimentally found until much later [2|. For BOs to be observable in an experiment, 
their periods have to be shorter than the scattering time r, and therefore the applied field 
has to surpass h/{elr). This value is too large for most natural materials, in which I is of 
angstrom size. In 1970, Esaki and Tsu suggested to create an artificial crystal, which they 
called a superlattice (SL), by growing many identical periods comprising a number of layers 
of two different semiconductors with similar lattice constants [3] . The period of the resulting 
one dimensional crystal may be much larger, say about 10 nm, and this gives reasonable 
electric fields of about 10 kV/cm, which are within the range of experimental observation. 
Damped Bloch oscillations were first observed in 1992 in semiconductor SLs whose initial 
state was prepared optically In recent years, BOs have been observed in other artificial 
crystals such as atoms placed in the potential minima of a laser-induced optical standing 

n n 

wave [4|, photons in a periodic array of waveguides [5| and Bose-Einstein condensates in 
optical lattices [§] among other systems 

BOs are potentially important to design infrared detectors, emitters or lasers which can 
be tuned in the Terahertz frequency range simply by varying the applied electric field 
Another application is based on the fact that BOs give rise to a resonance peak in the ab- 



sorption coefficient under dc+ac bias and a positive gain at THz frequencies 



The latter 



has been observed in quantum cascade laser structures {9J. These applications are severely 
limited due to scattering which rapidly damps BOs and, for a dc voltage biased SL, favors the 
formation of electric field domains (EFDs) whose dynamics yields self-sustained oscillations 
of lower frequency (GigaHertz) 10, ll| (a phenomenon similar to the Gunn effect in bulk 
GaAs [3]). EFD formation may also preclude THz gain in simple dc+ac driven SL which 
is typically calculated assuming spatially uniform solutions of drift-diffusion or Boltzmann 
type equations [7, [131 16]. This assumption has not been tested by solving space-dependent 
equations with appropriate boundary conditions or by experiments in semiconductor su- 
perlattices. An interesting idea for efficient terahertz harmonics generation is to excite 
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relaxation oscillations in the superlattice by incident radiation from a waveguide [171 ]. 

To understand the role of EFD formation in the observation of BOs or THz Bloch gain, our 
starting point should be a model in which BOs and EFDs are both possible solutions of the 
governing equations. One simple possibility is to consider single miniband SLs described by 
Boltzmann-Poisson transport equations. In 1971, Ktitorov, Simin and Sindalovskii (KSS) 
considered a one- dimensional Boltzmann transport equation with a collision model com- 
prising two terms: a simple relaxation to equilibrium and another term describing energy 
conserving impurity collisions [8J. Later Ignatov and Shashkin replaced relaxation to a 
Boltzmann type local equilibrium proportional to the instantaneous electron density in- 
stead of relaxation to global equilibrium 18) . More general relaxations to Fermi-Dirac local 
equilibria and self-consistent coupling of electric field and electron density via the Pois- 



son equation have been considered recently 19j-[2l|. The characteristic equations of these 



Boltzmann-Poisson models exhibit BOs as solutions and there is a hydro dynamic regime 
for large applied electric fields described by drift- diffusion equations 19[. However both the 
Boltzmann-Poisson and the drift- diffusion systems do not have BOs as solutions. Instead, 
they display self-sustained oscillations of the current through the SL due to periodic nu- 



cleation and motion of EFDs [19] , |20j, similar to the Gunn effect in bulk GaAs [12]. The 
reason of this shortcoming becomes clear when the equations for the moments of the distri- 
bution function are analyzed. It turns out that the current density and the energy density 
oscillate at the Bloch frequency during BOs but the electron density varies slowly. Thus a 
local equilibrium that depends only on the electron density (and not on the instantaneous 
value of the current and energy densities) cannot produce equations for these magnitudes 
in the hydrodynamic limit. The situation is similar to that found in gases which motivated 
the Bhatnagar-Gross-Krook (BGK) collision model 22] . If we want to derive hydrodynamic 
equations for the mass density, average velocity and temperature of a gas using relaxation to 
local equilibrium as a collision term, the local equilibrium distribution must depend on these 
magnitudes [22]. In this paper, we study a collision model similar to BGK's, i.e., relaxation 



to a local equilibrium depending on the electron, current and energy densities [23]. The most 
important property of the proposed model is that it allows the local equilibrium distribution 
to oscillate at the Bloch frequency, which is the crucial feature (missing in the KSS kinetic 
equation) if we want to derive a hydrodynamic regime that allows BOs. Since the scattering 
processes in a SL dissipate energy and momentum, our collision model includes two nonzero 
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restitution coefficients. This is similar to the case of low density granular flows in which 



inter-grain collisions preserve momentum 
model includes one restitution coefficient 



3ut dissipate energy and the corresponding BGK 



24j. In low density granular flows, inelastic colli- 



sions dissipate energy and, as a consequence, the granular gas is cooling down continuously 
(unless there is a continuous injection of energy, for example through the boundaries). Nev- 
ertheless using a Chapman- Enskog method, it is possible to derive hydrodynamic equations 
about a local equilibrium with a temperature that is continuously decreasing, the so-called 
homogeneous cooling state [241 ] . The hydrodynamic equations derived in 24j for the simple 
dissipative BGK model of low density granular flows have also been obtained for the inelastic 
Boltzmann equation (with an integral collision kernel) in a double limit of small Knudsen 



number and almost elastic collisions 



25]. 



In this paper we derive and solve numerically hydrodynamic equations containing BOs 
and EFDs among their solutions for a dc voltage biased SL. Bloch gain for a dc+ac driven 
SL will be studied elsewhere. Hydrodynamic equations are derived in a double limit: (i) the 
field-dependent term and the collision term are of the same order and dominate all others in 
the kinetic equation, and (ii) the collisions are almost elastic so that energy and momentum 
dissipation are of the same order as the spatial gradients in the balance equations. Extensions 
of classical kinetic theory methods based on assumption (i), such as the Chapman- Enskog 
technique, yield transport coefficients which become singular if the electric field becomes 
zero. Fixing this shortcoming requires matching to a multiple time scales expansion based 
on assumption (ii) and on a quasi-stationary solution of the equations for the first moments 
of the distribution function. Our techniques might be useful in other problems in kinetic 
theory having a similar structure. 

Once these difficulties are overcome, we can show that, in the appropriate limit, the 
electron current density and mean energy oscillate at the Bloch frequency, whereas the 
electron density, the electric field and the envelope of the BOs vary on a slower scale and 
are described by balance equations (hydrodynamic regime). Appropriate boundary and 
initial conditions include initiation of the BOs possibly by optical means [2]. Numerical 
solutions in the appropriate parameter range show that initial profiles for the field and 
;he BO amplitude evolve to stable spatially inhomogeneous profiles at room temperature 



261 ] . At low temperature (70 K), we have found that Bloch oscillations and Gunn-type 



oscillations due to EFD dynamics may coexist. Increasing lattice temperature produces large 
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diffusion coefficients in the electron current density (averaged over the BOs) as compared 
to the convective part thereof. This eliminates the Gunn-type oscillations. At low lattice 
temperature, the diffusion does not change that much, but convection dominates the average 
electron current density, thereby facilitating movable EFDs and Gunn-type oscillations. This 
novel finding of coexisting BOs of about 0.36 THz and 13.8 GHz Gunn-type oscillations 



is somewhat unexpected as it is usua! 



ly assumed that Gunn-type oscillations have to be 



eliminated to get BOs or THz gain 

The rest of the paper is as follows. In Section [Til we review the KSS-Poisson transport 
equations. These equations are written in nondimensional form in Section UTTI and equations 
for the two first moments of the distribution function (electron, current and energy densities) 
are derived from them. The moment equations do not form a closed set because they 
depend on the second moment. Assuming that the second moment is a function of the 
lower moments to be found later (the closure expression), we derive a reduced system of 
modulation equations for electron density, electric field and amplitude of the BOs by means 
of a nonlinear multiple scales method. Its analysis shows that BOs are always damped for 
the KSS-Poisson model. In Section \1V\ we present our dissipative BGK collision model 
whose local equilibrium distribution depends on electron, current and energy densities. The 
corresponding modulation equations may support self-sustained BOs as solutions. They 
contain closure functions of electron, current and energy densities that have to be calculated 
by a different method. The precise form of the closure expressions are found in Section 
mby etching the equations for electnc field, electron density and anrplitude of the BOs 
found in Section II III to the result of using the Chapman-Enskog method (CEM) [27[ on 
a modified kinetic equation for a distribution function that depends on the BO phase. 
This method yields equations with transport coefficients which are singular in the limit of 
vanishing electric field but they are compatible with the modulation equations of section 
II I II in an intermediate limit of sufficiently small fields. This compatibility yields the sought 
closure expressions. The results of numerical simulations of the modulation equations with 
appropriate boundary and initial conditions are presented in section IVT1 Section IVlII contains 
our conclusions. IA1 shows that Bloch oscillations are always damped for the dissipative BGK 
model with some particular local equilibrium distributions. [B]gives some technical details on 
the local Boltzmann equilibrium. In[Cl we derive a drift-diffusion system for the electric field 
in superlattices with strongly inelastic collisions by using a CEM similar to that described 
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II. THE KSS BOLTZMANN-POISSON MODEL 



We shall present our ideas in the very simple case of a n-doped semiconductor SL having 
only one populated miniband with the tight-binding dispersion relation: 

Al 



£(k) = — (1 — cos kl), v(k) = — -tt = — sin kl . 
K ' 2 K h v ; h dk 2h 



(1) 



Here A is the miniband width, I the SL period, H is the Planck constant and v (k) is the 
electron group velocity. Electron motion and the electric field are directed along the SL 
growth direction which we take as the x axis. In this case, the following modified KSS 
model describes electron motion including impurity collisions (which conserve energy but 
not momentum) and inelastic electron-phonon collisions 191 ]: 
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Q P (f) 



f w (k;n) dk = n, 



- Vp Af=—*\f{x,k,t)-f(x,-k,t)]. 



(2) 
(3) 

(4) 

(5) 
(6) 

(7) 
(8) 



Here /, n, No, e, k B , — e < 0, m*, /i and — F = —dW/dx are the one-particle distribu- 
tion function, the 2D electron density, the 2D doping density, the dielectric constant, the 
Boltzmann constant, the electron charge, the effective mass of the electron, the electro- 
chemical potential and the electric field, respectively. W is the electric potential. Note 
that the ID distribution functions have the same units as the 2D electron density n and 
that the electrochemical potential \i is a function of n obtained by solving (E])-©- The ID 
Fermi-Dirac local equilibrium fl6]) is obtained by integrating the 3D Fermi-Dirac distribution 
1/(1 + e^~ £( ^l( hBT ° ) ) with £(k) = £{k) + h 2 k]_/ '(2m*) over the transversal wave vector k±. 
Tq is the lattice temperature, u e and v p are collision frequencies which we take as given 



constants. The distribution function is periodic in k with period 2ir/l. A quantum version 
of the semiclassical equations ©-dEJ) was studied in [28 ] . 

The KSS-Poisson system ([2])-© goes beyond relaxation to equilibrium and linear response 
theory. The collision terms in ([2]) push the distribution function close to the local Fermi- 
Dirac equilibrium which depends on the instantaneous value of the electron density as 
indicated by (J7|). In the case of a finite SL biased at zero volts and having insulating contacts, 
we can show that the system evolves toward a global equilibrium ([6]) with F = 0, n = Np 
and the chemical potential corresponding to this doping density by finding a free energy 



functional and using it to prove the H theorem 



211 ] . If the SL has Ohmic contacts and is 



subject to an appropriate dc voltage, Gunn type self-sustained oscillations are possible and 



the free energy oscillates at the same frequency 



2l|. 



A. Characteristic equations, moment equations and Bloch oscillations 



The characteristic equations associated to (j2J) are 

dx A/ dk eF 

— = vik) = — sin kl, — = — — , 

dt K J 2h ' dt h 



(9) 



which, for constant F, have BO solutions x(t) = — cos( £ P (t — t )) . The influence of 
scattering can be seen from the equations for the moments of the distribution function. Since 
S(k) and f 1D are even in k and v(k) and Af are odd in k, the collision operators Q e (f) and 
Q p (f) satisfy the conditions: 

pit /I 



where 



r Q e , P (f)dk = 0, r £(k)Q p (f)dk = 0, 

J-tt/1 J-tt/1 



e 

2^ 



■k/I 
-■k/1 



■k/1 
—k/1 

A 



-£{k) 



Q e (f) dk = -v e n{E - E 



1D\ 



J n (x,t) = 
E(x,t) = 
E 1D (x,t) 



2tt 



■k/1 



k/1 



v(k) f(x, k, t) dk, 



I 



k/1 



A 



2nn(x,t) J_ n/l [ 2 

i r /l 

2nn(x,t) J^n 



£{k) 



f(x, k, t) dk, 
f 1D (x,k,t) dk, 



(10) 

(11) 
(12) 

(13) 
(14) 
(15) 



are electronic current and energy densities. Thus Q e (f) dissipates energy and momentum 
whereas Q p (f) dissipates momentum but not energy. For a finite SL with insulating contacts 
and zero voltage bias, these collision terms dissipate the electron energy and momentum 
until the electrons reach equilibrium at the lattice temperature T , n = Np, F = and zero 
current 21]. 



To obtain equations for n, J n and E, we multiply ([2]) by 1, v{k) and A/2 — £(k), respec- 
tively, integrate over k and simplify the results by means of (fT0]) - fTl2"]) . thereby obtaining 

e On dJ n 

Tm + ^ = > (16) 

dJ n eAH d eHnEF 

dE_lEdJ^_^l d_ lmh + FJJ = , E _ ElD) (lg) 
dt en dx 8hn dx n 

Here we have used and the Fourier coefficients fj of the periodic distribution function: 



oo 



f(x,k,t)= fi&t)e ijkl - (19) 

j=-oo 

Note that J n = —eA lmf 1 /(2h) and E = ARe/i/(2n). We can eliminate the electron 
density from ( TT6l) by using the Poisson equation ([3]) and integrating the result over x, thereby 
obtaining the following form of Ampere's law 

e^ + J n = At). (20) 

Here J(t) is the total current density. Note that f|T6l) - f|T8|) are a closed system of equations 
in the case of space independent moments. The dissipation terms in the right hand side of 
f|T7|) and fTTSl) ensure that a global equilibrium / = f 1D with n = N D , F = 0, J n = J = 
and E = E 1D is reached 



2l|. 



Note that space independent solutions of (|T6|) - (Tl8|) with u e = v v = (elastic collisions) 
have a constant electron density n, whereas J n and E satisfy the equation of a linear oscillator 
with the Bloch frequency wsioch = eFl/fv. 

dJ n eHF BE IF d 2 J n e 2 l 2 F 2 

-df-l^ nE = ^ + - Jn = ^^ 2 - + — Jn = - 



S 



III. MOMENT EQUATIONS AND DAMPED BOS FOR THE ALMOST ELASTIC 
KSS MODEL 

In this Section, we shall derive equations for the amplitude of the BOs in the limit of an 
almost elastic KSS transport equation. We shall use a quite general perturbation method 
that will be applied to other Boltzmann transport models later in this paper. Our results will 
show that BOs cannot be sustained within the KSS model and point out to the insufficiency 
thereof. 



A. Nondimensional KSS-Poisson and moment equations 

To study the KSS-Poisson transport equations and its associated moment equations, it 
is convenient to nondimensionalize them using the units indicated in Table [H They are 
appropriate for the hyperbolic limit S — > 0, in which the collision and Bloch frequencies 
are comparable and the corresponding terms dominate all others in fl2]). Let v be a typical 
collision frequency related to v e and v v . The field-dependent term in (T5]) has the order 
e[F]l[f]/h, whereas the collision terms are of order v\f\. Here [/] and [F] are typical scales of 
distribution function and field, respectively. Equations ([3]) and (@J with [A;] = 1/7 imply that 
[/] = [n] = Np. Collision and field dependent terms are of the same order for [F] = hu/ (el). 
From the Poisson equation we obtain: [x] = = ^jf^- The ratio from the convective 
term proportional to [v(k)] = Al/(2h) to the collision term of order v is a small dimensionless 
parameter 



6 - (21) 



This is also the ratio between the scattering time and the dielectric relaxation time and it 
plays the same role as the Knudsen number in the kinetic theory of gases. Defining now 
/ = f/Np, n = n/No, E = 2E / A, J n = J/[J n ], x = x/[x], . . . (where [y] are the units in 
Tabled]), we can rewrite all equations so far written in nondimensional form. Omitting the 
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hats over the variables, we find the following nondimensional versions of ©-(jHl) 



dF 

dx 



n — 1, 



n 



1 

2^ 



f(x, k, t) dk 



2tt 



f lu (x,k,t) dk, 



J V ' ; irh 2 N D 



/2a- A A 
1 + exp I - 



2k B T 2k B T 



cos k 



(22) 

(23) 
(24) 

(25) 



x k 



eN D A ehv 1 2eh 2 u e 2 NplA 
2h i^/Vo 2 e 2 N D lA 2eK 2 v 2 



f,n F £, E v(k) J n 

]\T hv A IA 

D H 2 2h 

10 10 cm -2 kV/cm meV 10 4 m/s 10 4 A/cm 2 nm 1/nm ps 

4.048 130 8 6.15 7.88 116 0.2 1.88 0.0053 



TABLE I: Hyperbolic scaling and nondimensionalization with v = 10 14 Hz. 

The moment equations (fT6"l) - (fT51) in nondimensional form are 

dn dJr, 



dt ~^ dx 
nEF = 5 

FJ n 

n 



5 



0. 



dJn j_ 1 9 ( d nj t 

-W + 2Tx {n - Reh) + 1 ^ 

BE E dJ n Id 



dt n dx 2n dx 



— —lmf 2 + le (E - E 



1D\ 



provided we define 7 e and 7 3 - through the relations 



— = (fye, 

V V 



6jj. 



(26) 
(27) 

(28) 
(29) 



We can rewrite (|2Tj) and (1251) in terms of fi = nE — iJ n , the first harmonic of the 
distribution function: 

oo 

f(x,k,t;5)= ]T ])(,: /:<))< (30) 



3=-oo 



as 



where /i is the complex conjugate of /i. The Ampere law (T2"Uj) is simply 

<9F 



ft 



+ J„ = J. 



(31) 



(32) 
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B. Amplitude of the Bloch oscillations 



The moment equations (126!) and (I3TI) for n = fo and fi are not closed because the higher 
moment A appears in them. In general, equations for moments /„,..., / will contain 
terms depending on f n+ \. Singular perturbation methods, such as the CEM [27] . produce a 
closure expression 

f 2 = g(n,F,f l] 6), (33) 
in which g can be written as a power series in 5. For the KSS-Poisson model, such an 



expression was derived in [19|. We will assume for the time being that g is a given known 
function and derive modulation equations for the slowly varying quantities n(x,t), F(x,t) 
and A(x, t). 

If we assume (as it is usually done in the method of multiple scales) that the moments and 
the field are functions of both a fast time scale r — t/5 (corresponding to a dimensional time 
unit l/u) and the slow time scale t, n = n(x, T,t;5), F = F(x, r, t; 5) and fi = fi(x, r, t; S), 
so that dn/dt in fl26|) becomes dn/dt + d^dn/dr and so on. Equations fl26l) . f l3~Tj) - fl32]) 
should be replaced by 

*" + (34) 



dr \ dt dx 

9F . ( T T dF\ 

<M J - J n-^r , (35) 



dr V 
d \ 



Id dh 



(36) 



Setting now 5 = 0, we find 

n = n(x,t), F = F(x,t), fi = A(x, t) e~ iFr , (37) 

in which n(x,t), F(x,t) and the envelope function A(x,t) do not depend on the fast time 
scale. (13^|) indicates that n varies slowly on the time scale t. Similarly and according to (|35|) . 
F is independent of r provided the total current density J(t) is of order 1. In practice, the 
size of J is set by J n and by the bias condition. Imposing a voltage bias condition between 
contacts at the ends of a SL with finitely many periods, J = 0(1) if we assume that this 
voltage is constant or it varies on the slow scale t. We shall not consider in this paper the 
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case of voltage bias varying on the fast time scale r (e.g., an ac voltage biased SL driven at a 
frequency of order 1/5), for which J = 0(1/5), and we have to modify the present analysis. 

The solution (157)) of ( I34p - ( 13 6 j) for 5 = exhibits BOs with frequency F. Before deriving 
modulation equations, it is useful to get first a quasi- stationary distribution function that 
solves (155")) and is independent of r: 



fl,i 



5 



+ 



F 2 + 5 2 lj% 

hi 



-y e nE w (5j j -iF) + 



F + i5-f e d 
2 dx 



[n - Reg s ) 



— iF d 

2 dx lmg ' S + ^ F ~ Slj ^ Re h s~( F + ^7e)Im h s 



(38) 



provided we have replaced h(x, t) = dfi/dt. We introduce the function h(x, t) because extra 
terms having this form appear in the moment equations when we use the CEM. The specific 
expressions for g$ and hs will be obtained by matching our results in this Section with those 
obtained by the CEM. See Section IVl and [Cl 



Remark 1. All terms in Equation (138)) have nonzero limits as F — > and this equation is 
the key step in the regularization of the results obtained in Section |V] using the CEM. 



Eq. (138)) is equivalent to 



J, 



n,S 



F 2 + 5 2 1]le 



le E lD nF+^-^-lmg s 
2 ox 



(39) 



F 2 + 5 2 l3le 



5 l3 I le E lD + ^^lmg s 



1 d 



2n dx 
5j 
n 



tReh s 



Im hi 



77 



(40) 



The subscript S in g$ and in hs stresses that these functions are calculated with r- 
independent n, F, J n< s and Es- Note that for F = 0(1), fi } s = 0(5/ F), whereas f\ t s = 0(1) 
if F <C 5. Thus the order of f\ t s depends on the order of magnitude of F and it is better 
to treat the compact expression (138)) as an 0(5) quantity. Without the x-derivatives and 
t-derivative in the functions g and h, the right hand sides of ( 1391) and (j40p correspond to the 
uniform stationary state 



J n ,su = nv d (F;n), v d 



5 le E lu F 
<5 2 7i7e + F 2 ' 



E. 



Su 



5 2 lel^ D 

5 2 ljle + F 2 ' 



(41) 
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The Boltzmann limit of ( 125]) . corresponds to approximating ln(l + x) ~ x in that expression 
and then calculating the chemical potential by means of (124ft . The resulting expression 
jiD ^ n7re* c ° sfc // (A)), where j3 = A/(2k B T ) and Iq(x) is the modified Bessel function of 



index zero 



29l | . produces a constant value E 1D = /1 (0n)/In(0n)- Inserting this in the drift 



velocity, (j4ip gives the well-known Ignatov-Shashkin formula [18] 



2v u T Al I^Pq) r eFl ju e + u p 



l + J 72 4hT e I ((3 ) Were V 

which we have written back in dimensional units. It reduces to the Esaki-Tsu drift velocity 
in the limit 0o — > oo (zero lattice temperature), in which the Bessel functions are absent. 

Remark 2. Comparing the Ignatov-Shashkin formula fj42|) to experimentally obtained 
current-voltage characteristic curves yields the numerical values of the collision frequen- 



cies u e and u p [30]. 



According to (l4"0]) . the mean energy E decreases as the field F increases, whereas the 
average energy (8) obtained by averaging (pQ), 

(S) = -— / £fdk = --E, (43) 



2vrn 2 

increases with the electric field, as one would have expected. Note that (j4"Tj) is an asymptot- 
ically stable stationary solution of the moment equations f|T6|) - f|T8|) provided we ignore the 
spatial dependence of n, J n and E. 

If we insert fi = fi t s(x,t) + $(x, t,r) in ( )36l) . we obtain the equation: 



(44) 



2 2 ^ 2i9a; 

Since F and n are still varying on the slow time scale t, it is appropriate to introduce the 
following nonlinear fast time scale instead of r: 

1 f* 

6 = - F(x,s)ds, (45) 
o Jo 

which yields 89/ dt = F/5, 89/dr = F. Note that, in dimensional units, the phase 9 equals 
the integral of the Bloch frequency eFl/h over dimensional time, and therefore the partial 
derivative of 9 over dimensional time equals the Bloch frequency. Thus 9 is the phase of the 
Bloch oscillations. 
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The fast and slow time scales 9 and t will be used to set up a method of nonlinear multiple 
scales below in order to find out the modulation equations on the slow time scale t. If we 
consider n, F and $ to be functions of x, 9 and t, Eqs. (I3~4"]) . 035]) and f l4"4"]) become 



-5 



^-#Im (A 5 + $) 
9t 2 



(46) 



2 2z ox 

The method of multiple scales is based on the expansions: 

i 

n(x,t;5) = ^5 m n (m) (#,2,t) + 0{5 2 ), 

m=0 
1 

F(x, t;8) = J2 S m F {m) (0, x, t) + 0(5 2 ), 

m=0 
1 

$(x,t;S) = ^5 m $ (m) (^,x,t) + 0{5 2 ), 



(47) 



(48) 
(49) 
(50) 



771=0 



and on assuming that n^ m \ F( m ) and $( m ) are 2/T-periodic functions of 9. Inserting 
(150"]) in ( 146]) - (1471) and (123]) . we obtain the following hierarchy of equations: 

a n (o) 



89 
dF<® 
dx 



0. 

n<°> - 1, 







i = o, 



F { 0) ^r^- = — + 77- Im (/i, s + $ (0) ) 



dx 



dt dx 



n 



(i) 



09 dt 2 



$(o) 



(51) 
(52) 
(53) 

(54) 
(55) 

(56) 



and so on. 

The solution of (1531) is 



= A(x,t) 



(57) 
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whereas (15T1) and (|52p indicate that n^ ' and do not depend on 9 31]. The solutions of 
(15^1) and (1561) are 27r-periodic functions of # only if the right hand sides of these equations 
do not contain secular terms proportional to 1 and e~ ie , respectively. This is the case if the 
integral of the right hand side of f )54|) and the integral of e ld times the right hand side of 
( 156]) over [— 7r, 7r] are both zero. These solvability conditions give: 

d n (0) Q 



dt 2 



A + ^ f_f9(n^,F^, fl , s + Ae-^0)^. (59) 



Instead of (l58j) . we can use the Ampere's law (132]) averaged over 6* with (J„) = — Im/^s 
given by (1391 : 



+ 



dt F(°) 2 + 5 2 7j7e 



^V^ ^ — ^-Im^ 



|^ {0) - Re^) - F< >Re/i s + S le lmh 



2 9a; 

1 (J)e, (60) 



where (J)e is the total current density averaged over one period of 9. Equations (|59|) . ( 160]) 
and ( 152]) (the Poisson equation) describe the Bloch oscillations. 

It is important to note that ( 158]) and (160]) are continuity and Ampere's equations averaged 
over the fast scale 9. The total current density depends on the bias condition. For a dc 
voltage biased SL of nondimensional length L = (N + l)l/[x], we have 



- / F(x,t)dx = 0, (61) 
L Jo 

where = eV/[hu e (N + 1)]) is a dimensionless average field proportional to the constant 
applied voltage V . Integrating the Ampere's equation fl32l and using d<f>/dt = 0, we obtain 

3 = I jf JndX = I / L[Jn ' S " Im (Ae ~ i6)] ^ (62) 

where we have used f\ ~ /^g + Ae - * 6 and J n — — Im/^5 = J n) g— lm(Ae~ ie ), where J n> s 
and 9 are given by (139]) and (HSJ), respectively. Eq. ( 160]) and the dc voltage bias condition 
yield 

(J)e = \f J n ,sdx^J- (J) e = ~\ \ Im i A ^ W ) dx. (63) 
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C. Insufficiency of the KSS-Poisson model 



The local equilibrium ( 125]) depends only on n, therefore g is a function of n and F but 
it does not depend on A and 6. Thus = 0, the second term in the right hand side of 
( 15 9 p is zero and therefore A(x, t) = A(x, 0) e - ^ 4 " 7 ^*/ 2 . The amplitude of the BOs decays 
exponentially fast to zero. This is consistent with the previous result that the hydrodynamic 



limit yields only a drift- diffusion equation for n and the electric field 



19]. 



IV. DISSIPATIVE BGK COLLISION MODEL 



We have shown that the KSS-Poisson model cannot sustain BOs because its local equi- 
librium function does not depend on J n and E and therefore does not depend on the Bloch 
phase 6 when fl57|) is used. Equation f )60|) becomes a drift-diffusion equation in this case. 
Similarly to the original BGK collision model [22( , we need a local equilibrium distribution 
that depends on n, J n and E in order to obtain a richer set of hydrodynamic equations. To 
account for thermal effects, we replace the following more general Fermi-Dirac distribution 



instead of / 



ID 



23 



26] 



32| : 



/ A*ccj T a ^ — 
in dimensional units, or 

f 1Da (hJ,u,jl) 

with 

P 



rrfhsTa 



In 



7T 



1 + exp 



/i a + hku a — £{k) 
k B T a 



m*A 



2nph 2 N D 

jjL a flUr 



In 1 



8 



jl+uk—j3+p cos k 



A 



(64) 



(65) 



(66) 



ksT a ksTfJ, 2ksT a 

in nondimensional units. In (IM|) . hu a k should be considered a periodic function of k with 

period 2n/l. Then f 1Da is 27r//-periodic in k, same as the electron distribution function 
/. The multipliers fi a , u a and T a should be selected so that the electron density (jl]), the 
electronic current density ([TBI and the mean energy (T14"]) satisfy the equations: 

l_ Mi 
2n 



f 1Da dk 



n, 



e 

2^ 



TV/I 

■k/1 



—k/1 

I Mi 



2nn 



v(k) f iUa dk=(l-a 3 )J n , 

f 1Da dk = a e E + (1- a e )E. 



TV/1 



(67) 
(68) 
(69) 
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Here acj and a e are dimensionless restitution coefficients taking values on the interval [0, 1] 
(see below). Eq is the mean energy at the lattice temperature of the global equilibrium 
reached by a finite SL with insulating contacts and zero voltage bias. The nondimensional 
versions of (J6"7]) - (j6l?]) are 



— / f 1Da dk 



n. 



1 

— / sin k f 1Da dk 



1 



27m 



cos k f 1Da dk 



y 1 ^ j ) ^ n 1 

-- a e E + (1 - a e )E. 



(70) 
(71) 
(72) 



The nondimensional multipliers jl, u and are functions of n, J n and E determined by 
solving ( 17D|) - (i72l . With these definitions, the collision operator Q e (f) of (jSJ) with f 1Da 
instead of f 1D satisfies 



w/l 



w/l 



Qe(f) dk = 0, 



— / v(k)Q e (f) dk = -v e ajJ n 
2tt J..-H 



2im 



n/l 
w/l 

■k/1 



-S(k) 



Q e (f)dk = -v e a e (E-E ) 



(73) 
(74) 
(75) 



In nondimensional units, we have 



lDa\ 



dk = —a e n(E — Eq) — ajiJ n , 



(76) 



instead of (j74[) and fj75|) . provided we select v = u e as our unit of collision frequency. The 
restitution coefficients a.j and a e , < aj^ e < 1, measure the fraction of momentum and of 
energy lost in inelastic collisions, and correspond to the single restitution coefficient used in 



granular gases 



24j. Obviously for a e j = the collisions are elastic. Note that we do not use 



the temperature T a = aT as in granular gases because the relation between energy density 
and temperature is not linear in the present case. To simplify matters, we shall assume 
that the restitution coefficients are constant. For space independent solutions of the kinetic 
equation, this leads to exponentially fast decay of the average energy and momentum in 
contrast with the algebraic decay of energy found in granular gases 24]. 
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A. Choice of local distribution function 



The distribution function fl64|) has the same form as the equilibrium distribution for an 
electron-phonon collision term in which the phonon distribution is n q = (e^ R ( a, « -q ' u ) — 
where u q is the phonon frequency corresponding to a wave vector q and u is the average 



velocity. In fact, the electron-phonon collision term for a bulk semiconductor is [33 1 
C[f}(k) = [ K(\k - k'|) {5(8' -£ + hw q -hu-q) [n q f'(l - f) 

J IB 

-(1 + n q ) f(l - /')] + 5(8' - E - hw q + hq • u) 
x [(1 + „,) f'(l - /) - n q f (1 - /')]} dk>, (77) 
n * = effliK-u-q) _ i » q = k-k'. (78) 

Here depends on the phonon type, B is the Brillouin zone and / = /(k), /' = /(k') 

with similar notation for the dispersion relation 8(k). For a kinetic equation c^/ + v(k) • 
V x / = if s is a function of k and / with < / < 1, we find 



Of 



[ s(k,f)dk + V x - [ v(k)s(k,/)dk=- / / KT(|k-k' 
Jb jb 7b Jm 



I) 



X 



<$(£' -S + hLU.-hq-u) e -^-" ku )(l + - /) (1 - /') 



1-f 1-f J \df df 

The right hand side of this equation is always less or equal than zero for 
ds { ePP-te-") f\ 

— = In I — ) , i.e. for an entropy density, (80) 

df \ 1-f J 

s(k, /) = / In / + (1 - /) ln(l - /) + 0(8 - Hk ■ u)/, (81) 



and the corresponding integral over A; is a Lyapunov functional for homogeneous distribu- 
tions. (Note that s > (38 f — 1, and the corresponding integral over k is bounded below 
because both the energy and the volume of the Brillouin zone are finite). This is the H- 
theorem for the electron-phonon kinetic equation. From (1821) . we see that the corresponding 
equilibrium solution satisfies e^ £ ~ ftk ' u )//(l — /) = (independent of k), and therefore 

•W 4 ) = i + e /3[£(k)-/ik-u-/ t ] ' ( 82 ) 

is the equilibrium distribution. If we use (182|) for a SL, then the dispersion relation is 
8(k) = 8(k) + h 2 k 2 ± /(2m*), and integration of fl82|) over the lateral wave vector kj_ yields 
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a distribution function of the same form as (|64|) provided u is directed in the SL growth 
direction. 

A different choice of local distribution function could be 



nh 2 



111 



1 + exp 



/i Q £(fc k a ) 
k B T a 



(83) 



instead of (164"|) . The multipliers would now be T a and fc Q (instead of u a ), to be selected 
so that conditions (|671)-(l69l) be satisfied. The choice of k a would correspond to the wave 
packet momentum in Lei's formulation 34(. For the tight binding dispersion relation ([Q), 
substitution of cos(A; — k a ) = cos k cos k a + sin k sin k a in ( 183|) would yield a distribution of 
the following type: 



f (^) i^ai Pai T a 



rffksTa 
irh 2 



In 



1 + exp 



fMa - E{k) - v{k)P a 



(84) 



The distribution (I8^|) could have been obtained from the maximum entropy principle as 
suggested in [35| and, for the tight binding dispersion relation, it becomes (1831) selecting 
appropriately the multipliers fi a and T a ; cf. |A] We show in |A] that the BO solutions corre- 
sponding to transport equations that have local equilibrium distribution functions ( 1831) and 
(with tight binding dispersion relation) are always damped. 



B. Equations of the model 

Since (|74l) and ( 1751) show that our collision model dissipates both momentum and energy, 
we propose a simpler equation for the distribution function with Q p (f) = in ()2]) and ( 164)) 
as the local distribution function instead of (Q. Recapitulating, the equations governing our 
inelastic BGK model are ^ and © with Q p = and Q e = -v e {f - f lDa ) given by (J64J) 
and ( 1BT|) - (16T?]) . n, J n and E are the moments of the distribution function given by (jlj), ( TT3|) 
and ( 1T4|) . respectively. In nondimensional units, the equations of the model are ( 122|) - (f24"|) 
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with Up = 0, v = v e and f 1Da given by (I65I) instead of f 1 . 



F dk +f ~ f 

OF 

dx 



■ ,'df . , 9/ 
-d I — + sin k — 

at ox 



n — 1, 



n = 
1 

27 



1 

2^ 



f lUa {k ]n ,E,J, 



7T 

e ik {f ~ f 1Da )dk = 
m*A 



1 

2^ 



f 1Da (x,k,t) dk, 
-5j e n(E - £ ) - ^7i^n, 

111 | 1 + e A+ sfe -/3+/3 cos fe 



(85) 

(86) 
(87) 

(88) 

(89) 



2irfih 2 N D 

Here we have assumed that the restitution coefficients are of order S and defined 7 e and jj 
for this model as: 

a e = 57 e , a j = djj, (90) 

instead of (129]) for the KSS equation. The restitution coefficients a e j can be fitted to 
experimentally obtained current- volt age characteristic curves in the same way as the KSS 
collision frequencies v e and v v . We find va e = v e and uatj = v e + z/ p . When modeling a finite 
SL, we need boundary conditions for / and F at the contacts attached to the SL boundaries 
and an initial condition for /. See References 11[ and 20J for a discussion. 



C. Boltzmann distribution 

We can simplify the previous formulas in the low temperature limit in which (3 — > 00, 
u = 0(j3), ft — > — 00 in (I89p . which becomes the Boltzmann distribution 

27rH 2 pN D 

and integrals over are calculated using Laplace's method. For sufficiently high temperature, 
the Boltzmann distribution fl9Tl) is again a good approximation and it yields simpler formulas. 
The parameter ft can be explicitly calculated using (I9"T|) in ( 157)) and the resulting distribution 
is 

„uk+B cos fc 

fB = „ ?1£ (92} 

eP<*»Kco8k(uK)dK' 
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in which u and $ are obtained in terms of J n /n and E by solving (IBB1) . As shown in[Bl the 
latter equations yield 

(1 — ctj) J n u e _/3 sinh(w7r) 



n 



P P Jo ^ c ° s K cosh(uK) dK 



P 



(93) 



E - a e {E - Eq) 



Jo 1 



? cos K 



cos K cosh(-ufT) cLfT 



p 8 cos K 
C 



cosh(-uK) rfK 



I (P) + 2u*YZx tw 1 ^) 



(94) 



where I s (x) are modified Bessel functions (29[]. At the lattice temperature, (3$ = A/(2A;sTo 
and for zero current, tt = 0, E = E , and (|94|) yields 

h0o) 



Eq 



Wo) 



(95) 



Further simplification follows if we impose aj = 1 in (1931) (which implies u = 0) so that the 
BGK collision term dissipates momentum and energy according to (|88|) . Then (|92|) becomes 

3/ 8 cos fc 



/oC9) 



(96) 



and /3 is obtained in terms of E by solving (1881) with u = 0, i.e. 

hiP) 



a e E + (1 - a e )E 



(97) 



The Fourier coefficients of the Boltzmann distribution ( 196]) are simply 



f B 

J 3 



2tt 



e~ ijk f B {k-n) dk = n 



(98) 



V. CHAPMAN-ENSKOG METHOD FOR ALMOST ELASTIC COLLISIONS 



In this Section, we shall derive modulation equations for n, F and A in the case of almost 
elastic collisions with < a e j 1, more precisely in the double limit of small "Knudsen" 
number 5 and vanishing restitution coefficients a e j = Sjej. For the case of granular gases, 
Sela and Goldhirsch derived hydrodynamic equations from the inelastic Boltzmann equation 
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using a CEM in a similar double limit 25] . We shall use the CEM [27| to obtain modulation 
equations for the electric field, the electron density and the amplitude A of the BOs. Then 
we will compare these equations with (I59p and (loTij) and identify g and h. 

We can repeat the calculations of Section HTT1 with the local equilibrium distribution ( I89p 
and get the same modulation equations (1591 - (ISOl) except that E replaces E 1D in them and 
7 e j are defined by ( 190]) instead of ( )29|) . Now the local equilibrium distribution f 1Da depends 
on fx = nE — iJ n and therefore it depends on the Bloch phase 9 through (J57J). Then the 
second term in the RHS of ( |59|) is no longer zero and BOs are not necessarily damped by 
scattering. This term will be now identified by using the Chapman-Enskog method to derive 
modulation equations that will be matched to (158j) - (loT?j) . To implement the CEM, we assume 
that the distribution function / is a function of 9, k, x and t, which is 27r-periodic in 9 and 



in k and that F is of order 1. Eq. (185]) becomes 



= + sm,g), (99) 

Cu(k, 9) = f(^- + ^-] u(k, 9) + u(k, 9), (100) 



89 dk / 

Equations (I99l) - (ll00p with F = 0(1) display a dominant balance between the collisions, 
the force due to the electric field and the change of / over the fast time scale 9. We are 
ignoring a possible fast relaxation from an initial distribution function to the BO distribution 
f(x,k,9,t;S) that depends only on one fast time scale 9. We are imposing the condition 
that / be periodic in 9 and considering the possibility of slow modulations of the BOs in 
the time scale t. 

The moment f\ = j\s + $, $ = Ae~ ld + 0(5), has a dominant part of order one, Ae~ ie , 
and a remainder of order 5. The remainder vanishes as 5 — > and it can be chosen not to 
contain a term proportional to e~ %e . Thus we assume: 

f l = Ae- W + 5B + 5 2 C + 0{5 3 ). (101) 

The local equilibrium f 1Da is a function of k, n and f\ through (1H71) - ( 1891) . Due to f llOll) . 
f lDa is a 27r-periodic function of k and of 9, which also depends on the slowly-varying 
functions n(x,t), F(x,t), A(x,t), B(x,t) and C(x,t). 
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Using the Fourier series 

f D «(k,9;6) = Y,fl? a e lUk+W \ (102) 



f{k,e ]X ,t,S) = Y,fjMf,s)e ijk+m , hi = r fe- ijk - iW f^%, (103) 

in (199]) with 5 = 0, we immediately obtain 

flDa(O) 

where the superscripts (0) refer to having set 5 = in ( 199|) (see below). 

The CEM starts from a leading order expression for the distribution function, (I104p . 
which does not depend explicitly on x and t. Instead, it depends on k and 6, and it is a 
function of quantities that vary slowly with x and t (n, F, A, B, C, . . .). B, C, . . . are to be 
calculated in terms of A, n, F and their differentials. While A, n and F are not expanded 
in powers of 5, their partial derivatives with respect to time (and therefore their equations 
of motion) are expanded instead. Thus the Chapman-Enskog Ansatz is 

oo 

f(x, k, t;S) = J2 f (m) 0; F, n, A B, C) 5 m , (105) 

m=0 

8F °° 

— + J im) (F n, A, B, C) 5 m = (J) 9 , (106) 

m=0 

^ = -J2^ {m) (Fn,A,B,C)5 m , (107) 

m=0 

— = J2^ m \F,n,A,B,C)5 m . (108) 

m=0 

We have used the Poisson equation (18"6"1) and (11061) to obtain (11071) . The local distribution 
function f 1Da can be expanded in powers of S, 



jlDa _ jlDa(m)fim^ (109) 

m=0 

and then flH7|) - flHBl yield the following compatibility conditions: 

f^ ] = flT {m) = n ^ (no) 

ft\ = Ci im) = A5 o m , (HI) 
fS = B, f^ a{1) = B + le nE , (112) 
/$ = C, fl Da{2) =C- le n ReB - i 7j lmB, (113) 
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and so on. Inserting ( 1 1 5 [) - (1 1 9 [) into (|99|) . we obtain a hierarchy of linear equations for the 
f( m ) whose right hand sides contain the functionals and A^ m \ The latter are calculated 
in such a way that the compatibility conditions flllOp - (1113)) hold. 
The equations for /W and are: 



£j(l) _ jlDo(l) 
£f(2) _ jlOa(2) 



sinfc4-), (114) 



dt 



+ sin k 



,, cte / dt 



(115) 



The subscripts m — 0, 1 in the right hand side of these equations mean that dF/dt, dn/dt 
and dA/dt are replaced by ((J)e5 m o — J^), —dj^ m '/dx and — -4.( m \ respectively. 
Upon insertion of (1104)) in (1114)) . the compatibility conditions ( II 10)) - (11131) yield 

(116) 
(117) 




+ J_^_( n _Ji^ ). (|| S) 

Note that B becomes singular in the limit as F — > 0. This is not surprising: we have 
assumed in this Section that F = 0(1) as 5 — > so that 6^0 and the Fourier series in 9 
of fi in (1101)) has a first harmonic Ae~ %e that is different from all other ones contained in 
B, C, etc. If F tends to 0, then the first two terms in (1100)) are smaller than the third one 
and the assumption (1 1 1 j) breaks down. Despite this shortcoming, we shall use the CEM to 
identify the closure functions g and h introduced in the previous Section. The coefficients 
in the resulting modulation equations are no longer singular. 

The compatibility conditions (11131) for provide the following functionals 

J (1) = -Im B, (119) 

^» = i%i. (120) 

Then the averaged Ampere's law and the equation for A including up to 0(5) terms are 

^--5\mB = (J) g , (121) 
£ = > + 7 ,M + i (122) 
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in which B is given by f II 1 8 j) and the Fourier coefficients of the solution of (11141) are 



./ 



(i) 

3,1 
(1) 



r (l) 
0,1 



1 + iF(j + /) 

jlDa(l) i ^ 



Of 



+ sink 4~) /(°). 
o dx J 



(123) 
(124) 



A. Identification of g and hs in the modulation equations 

We now impose that Equations ffmj) and (1T22]) match ([60]) (with £ instead of E 1D ) 
and ( |59|) . respectively. The result is that these equations match term by term in the overlap 
region 

1, 



(with B5 ~ fx t s) provided 



S'S — /2,o > 



(125) 
(126) 



Both equations hold if 



S = / 2 (0) +^ 



(i) 



(127) 



We have not yet calculated hs in ( J60l . To determine it, we require that the resulting 
equation for the field coincide with the drift- diffusion equation ( 1C.36I) derived in O for the 
case of inelastic collisions (without BOs). As seen in [Cj hs in (1601) should be replaced by 
the uniform part of dfi^s/dt in (1381) . i.e., 



d_ ( 5i e nE (6'Y j - iF) 
dt V 5 2 j e jj + F 2 



5j e E (5"/j - iF) dJ n ,Su 



8 2 lel] + F 2 



dx 



/ /n t x 5 SlenEoiSij - iF) 

+ {{J)g - Jn t Su) "rTTT 



<9F V ^ 2 7e7j + i 72 
Here J ntSu = 5^ e nE F '/(8 2 jjj e + -F 2 ) according to ( I4~T1) . 



128) 
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B. Modulation equations 



After straightforward calculations, we obtain the following reduced equations: 



dF 5 



dt F 2 + 5 2 ljle 



TP TP i F 9 T /2,0 

•y e E nF + — — Im 



lDa(O) 



tl("- Re rS5i- FRefts+ ^ Im ^ 



2 <9x l + 2iF 

: (J)o, (129) 



dt~ 2 A+ 2K&^ 1 + J' (1 ' )0) 

d dJ„ Su d \ f^-x 



. I. _ f l£»a(l) / .(0) 9 , // A t x 5 5 \ 



in addition to (I128p and to the Poisson equation (186]) : 

dF 



gx n-l. (132) 



To calculate /^-i^ in (11221) . we need to use 

flDa(l) _ f~(l) 9 , ~(1) g , 5(1) g ^ f lDa(0) fiQQ) 

In this expression, we should substitute jl^ \ u^°> and (3^ given by simultaneously solving 

/ 1Da(0) =n, /f^^ie-* (134) 

and also the solutions u^ 1 ', $^ of 

fl Da{l) = 0, (135) 
fl DaW = le nE + / 1)S - A e"* - Ze* e . (136) 

When we substitute these solutions, ( 1133}) becomes a function of fc, 9, n, A and /^g ~ B5 
which is 27r-periodic in k and 9. Its Fourier coefficient /2 -i * s then inserted in (11311) . 
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VI. NUMERICAL SOLUTIONS OF THE MODULATION EQUATIONS FOR n, F 
AND A 



A. Modulation equations 

According to the results in the previous section, the modulation equations for the local 
equilibrium distribution f lDa are (I129p . f 11 3 j) and HI 32j) . The coefficient functions appearing 
in these equations are given by fl!28j) . (113 ip and fl 133ft - ( TT36]) . The modulation equations with 
r^L-y = and lattice temperature of 300 K were shown in Ref. 



261 ] to exhibit BOs confined 



to a portion of the SL provided (j e + 7j)/2 < 7 C (j c is a critical value). 



B. Boundary, bias and initial conditions 

The boundary and bias conditions are: 

dF dF 
(J)e--^ = °vF (atar = 0), (J) e - — = a 1 nF (at x = L), (137) 



A = 0, at x = and at x = L, (138) 
F(x,t) dx = <f>. (139) 



1 ' L 



Lj 

The last equation is the nondimensional dc voltage bias condition ( IBTj) with = eV/[hu e (N+ 
1)]. It is also possible to set dA/dx = in the contacts and the numerical results are similar. 
According to (l6"2l and (1531 . the total current density is 

1 f L 1 f L 

J = T Jndx = T \- J ^ s ~ Im ( Ae ~ te )] dx 
L Jo L Jo 

i r L 

= (J)e - T / Im (Ae- W ) dx, (140) 
L Jo 

for dc voltage bias. In our numerical solutions, we have adopted uniform profiles for F(x, 0) 
and A(x, 0) as initial conditions. 

C. Numerical results 

We shall illustrate our results with numerical solutions of ( I129p - (ll32p with the Boltzmann 
local distribution function (|92|) -( 194|) . Using the more general Fermi-Dirac local equilibrium 
fl65|) complicates the numerical procedure by having to calculate one more multiplier at each 
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time step: fi in addition to (3 and u. We shall use indicative values similar to those in Ref. 30]: 
I = 5.06 nm, A = 16 meV, v P = 10 14 Hz, with a P = a, = 0.006 so that v P a P = v P ct; = 6 x 10 11 



Hz. 



tj = u.uub so tnat v e a e = u e aj 
The 3D doping density is N 3D = 8 x 10 16 cm" 3 , so that N D = N 3D l = 4.048 x 10 10 



cm -2 , and e = 12.85 e . We find 5 m 0.0053 and j e> j = a e j/5 = 1.1269. We consider 
a 50-period (N = 50) dc voltage biased SL with lattice temperature 70 K. We have used 
contact conductivities o"o = 60.6 (f2m) _1 and o\ = 15.15 (fim)^ 1 which yield dimensionless 
conductivities do = 1 and o\ = 0.25 in (11371) (conductivity units are [a] = e 2 Nr ) Al/(2h 2 i' P )). 
Initially, the profiles of A and F are uniform, with a common value 0.5501. 

For V = 0.133 V (therefore <fi = 0.04) and after a short transient that depends on the 
initial conditions, we observe coexisting BOs with frequency about 0.4 THz and Gunn type 
oscillations with frequency about 14 GHz. See the movie in the Supplementary matterial 
37|] . Fig. [1] shows several snapshots of the field and \A\ profiles of the Gunn type oscillation. 
While the amplitude of Gunn-type current oscillation is about 0.03 in nondimensional units 
(as seen in Fig. QJa) for the total current density averaged over the BOs), the BO part of 
the current oscillation has a larger amplitude of about 0.2; see Fig. |2]Ja). Fig. [2] illustrates 
the total current density (11401) of the coexisting 360 GHz Bloch and 13.8 GHz Gunn type 
oscillations, respectively. For each lattice temperature, there is a critical curve in the plane 
of restitution coefficients such that, for (j e + -fj)/2 > r y rr \t, BOs disappear after a relaxation 
time but they persist for smaller values of [pf e + 7j) 26]. 

Figure [3] shows the profiles of F and A and Fig. @] depicts the total current density at 
temperature 300K for the same values of a e j and the other parameters. We find BOs but 
not the slower Gunn type oscillations. Whether Bloch and Gunn type oscillations coexist 
depends on the relative size of the diffusion and convection terms in (11291) and ( 11301) which, 
in turn, are controlled by the lattice temperature. If diffusion terms are sufficiently small 
compared to convective terms in (11291) and (11301) (which happens for small enough lattice 
temperature), Gunn type oscillations mediated by EFDs in the F and A profiles are possible. 
For larger temperatures, Bloch and Gunn type oscillations cannot occur simultaneously. 
This latter fact was previously revealed by solving numerically a simpler version of the 
hydrodynamic equations with = in (11301) and quite large voltage bias |26j ]. Note that 
the largest peak in the current spectrum occurs at a lower frequency (260 GHz) than in the 
case of lower lattice temperature of Fig. 12(b). 
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FIG. 1: (a) 0-averaged total current density vs time during coexisting Bloch and Gunn type 
oscillations at 70 K. (b) Field profile vs space at the times t% to t$ marked in (a), (c) Same for the 
complex BO amplitude profile. To transform the magnitudes in this figure to dimensional units, 
use Table U [A] = N D . 
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FIG. 2: (a) Total current density vs time during coexisting Bloch and Gunn type oscillations 
at 70 K. (b) Fourier transform of the total current density showing two peaks corresponding to 
coexisting Bloch (0.36 THz) and Gunn type (13.8 GHz) oscillations. The zero-frequency constant 
corresponding to the time average of the total current density has been subtracted. 



VII. CONCLUDING REMARKS 



We have proposed a Boltzmann-BGK kinetic equation for electron transport in miniband 
semiconductor superlattices. Its local equilibrium depends on electron density, mean energy 
and current density and therefore it oscillates periodically in time with the Bloch frequency 
when the mean energy and the current density do the same. This model is richer than 
the usual BGK models traditionally used in this field and its corresponding hydrodynamic 
equations may exhibit Bloch oscillations which are absent in the hydrodynamic regime of 
the KSS and related models. We have introduced novel singular perturbation methods to 
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FIG. 3: Modulus of the BO complex amplitude and field profiles vs space for the stationary state 
at 300K. 



derive hydrodynamic equations describing Bloch oscillations in the limit in which collision 
and Bloch frequencies dominate all other terms in the kinetic equation and the collisions 
are almost elastic. By numerically solving the hydrodynamic equations with appropriate 
initial and boundary conditions, we find that nonlinearities may stabilize Bloch oscillations 
if the restitution coefficients are small enough. There are different scenarios depending on 
the lattice temperature. For sufficiently low temperature, Bloch and Gunn type oscillations 
mediated by electric field, current and mean energy domains may exist simultaneously for 
appropriate voltage ranges. These oscillations are spatially inhomogeneous and have field 
profiles with EFDs typical of Gunn oscillations. For larger lattice temperatures, Bloch and 
Gunn type oscillations do not coexist: the profiles of the electric field and the amplitude 
of the Bloch oscillations are independent of time but inhomogeneous in space 26[]. As the 
collisions become more inelastic, the parameter range for which BOs appear shrinks and 
these oscillations disappear for the standard superlattices used in experiments [30]. In the 
absence of BOs, the hydrodynamic equations become the known drift-diffusion system valid 
for inelastic collisions that may exhibit Gunn-type self-sustained oscillations due to periodic 
recycling of charge dipole domains for appropriate parameter values [11 ]. 



31 




100 200 300 400 500 
Frequency (GHz) 



FIG. 4: (a) Total current density vs time during Bloch oscillations at 300K. (b) Fourier trans- 
form of the total current density showing only one peak corresponding to BOs (0.27 THz). The 
zero-frequency constant corresponding to the time average of the total current density has been 
subtracted. 
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Appendix A: Local equilibrium distributions that cannot sustain Bloch oscillations 



Let us consider the local equilibrium distribution ( 183)) which becomes 

m*A 



f lD «(k;n,E,J n ) = m * In (l + eA-^-C*-*-)), ( A .l) 



2tt/3H 2 N d 

written in nondimensional units. Inserting this equation in ( 187)) . we find 



n = f} Da = m * A \n(l + e^ + ? cosk )dk, (A.2) 

27rf3h 2 N D V J 

after shifting the integration variable k — )■ (A; — A; a ). Similarly, we find 
3 2nBh 2 N D J-* V / 



m*A 



/7T 
cos(jA;) ln(l + e ^ + ^ cos/c ) dfc. (A.3) 



As 5 — >• 0, the left hand side of (1A.3)) for j = I becomes Ae~ ld according to (11121) . from 



which we obtain 



_m* /■ ln A +e A-/3+/3cosA cosfcdfc (A. 4) 

fc a = 0-arg(A). (A.5) 



Equations ( 1A.3j) and ( 1A.4I) can be solved to produce the leading order approximations of \x 



and (3 as functions of the slowly varying quantities n and \A\, whereas (1A.5I) indicates that 
k a varies rapidly as a shifted Bloch phase. Then (1A.3I) implies that fjf a ^ for I = —j and 
all the other harmonics are zero. We find that f^ {0) = = in and therefore 

BOs are always damped for this model. 

Let us consider now ( 1841) which, for the tight binding dispersion relation, can be rewritten 

as 

f 1Da (k;n,E,J n ) = m * A In (l + e i>*-0*+P*™ s k-P aS ink)\ (A 6) 

Replacing k = k a + £, we obtain 



(3* cos k — P a sin k = ((3* cos k a — P a sin fc Q ) cos £, (A. 7) 

provided 

P Q 

/3* sin /c a + Pq, cos k a = ==>- tan fc a = — (A. 8) 

(3* 
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Substituting HA.8|) in f!A.7j) . we get 



P 2 

f3* cos k — P a sin k = ( (3* + J cos k a )cos £, 



which, inserted in flA.6j) yields (1 A. 1 1) with 



[3* 



This shows that ([BID i s equivalent to 



(A.9) 



(A.10) 



Appendix B: Boltzmann local equilibrium distribution 



In nondimensional units, the Boltzmann distribution f[9"Tj) satisfying /q = n is 

uk+/3 cos fe 



7re" 



(B.l) 



fJV^COShfuit) 

The first moments of this distribution can be used to calculate {3 and u in terms of E and 
J n by solving 

g T e^ coaK cosh(uK) cos K dK 



J*eP cosK cosh(uK) dK 
£ cos K smh{uK) sin K dK 



a e E Q + (1 - ct e )£, 



/ 'V cosA ' cosh (uAT) dif 
The left hand side of (IB.3j) can be simplified by integrating the numerator by parts: 



(B.2) 
(B.3) 



e " sinh(ii7r) 



13 (3 eP cosK cosh(uK)dK 



(B.4) 



Equations (IB .21) and (1B.4|) contain the integral e^ cosfc cosh(uk) dk which can be calculated 
using the generating function [29] 



i cos k 



W) + 2j2li(P)cos(lk), 



(B.5) 



i=i 



with the result 



e^ cosfe cosh(w£;) dk 



u 



00 f-lV 



sinh(-u7r). 



(B.6) 
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From this formula we obtain 

-X In / e^ cosfe cosh (uk)dk = — ^hlZ±Ei 1^ 1 

We now use (1B.6j) and (IB. 7ft in (1B.2[) and (IB. 41) . thereby obtaining 



A(^) + ^ 2 E^fc^[^i(/3) + /m(/5)] 

— = a e £ + (1 - a e )E, (B.8) 



u 



/oC9) + 2fi=»E"i fcgW) 



'l-a 7 -)— . (B.9) 
n 



Appendix C: Inelastic collisions and the hyperbolic limit 

Here we shall use the CEM to obtain equations for the electric field and the electron 
density in the case of inelastic collisions with < a e j < 1. In the method of multiple scales, 
we expand the distribution function and all its moments in powers of 5 and consider slow and 
fast time scales. The condition that the terms in the distribution function be periodic (or, 
more generally, bounded as the fast time tends to infinity) in the fast time determines the 
modulation equations in the slow time scale. In the inelastic case, the damping coefficient 
(a e + 0£j)/2 in the equation for the BO amplitude is of order one. Thus the distribution 
function relaxes exponentially fast to a quasi-stationary function whose current and energy 
densities are given (to leading order) by (14T]) . This distribution is the starting point of the 



CEM which, in the inelastic case, is similar to that described in |19| and [27]. 

The leading order expression for the distribution function depends on time only through 
the moments n and F which vary on the slow time scale t. These moments are not expanded 
in powers of 5. Instead, their evolution equations are expanded (as we show below), and the 
corresponding terms in the expansion are determined so as to keep compatibility conditions 
issuing from the assumptions for the distribution function. The CEM can be used to obtain 
reduced equations for the moments containing terms of different order in 5, and this is 
something that the method of multiple scales cannot deliver. 

The leading-order distribution function is the solution of Eq. ( )82|) for 5 = 0. Its Fourier 
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coefficients are 

Re/f = h 1 + J fF2 h , (C.l) 

(0) Im fj Da — jF Re f} Da , x 

Wf = ^ 1 + ; 2F2 Jj ■ (C.2) 

We assume that J n and have already acquired their quasi-stationary values after a fast 
decay on the time scale r. These quasi- stationary values are functions of n, F and 5 to be 
determined now. The Chapman-Enskog Ansatz is 

oo 

f(x,k,t,6) = Y t f (m) (k;F,n)S m , (C.3) 

m=0 

^. + ^ |7 (m)( F>n ) ( J TO = J( t ), (C.4) 
m=0 

t = -^T x j{m) ^ n)5m - (C - 5) 

m=0 

In (jC.4p . the total current density is of course the same as its average over one period of 
the BOs. We have used the Poisson equation fl86|) to obtain (lC.5j) . The local distribution 
function f 1Da is now a function of n and F because J n and E depend now on n, F and S. 
We have 

OO 

j?lDa ^ ^ j?lDa(m) 

m=0 

and then (157)1 - (1551) yield the following compatibility conditions: 

/o (o) = f io am = (C7) 

ReA (0) = nE<°>, R-e/ 1 1 ' Da(0) = n [a e E + (1 - a e )£( )], (C.8) 
Imjf^-jf, ImA 1Da(0) = -(l-a,)4 0) - (G.9) 



Let us now find f 1Da ( m ) in (|C.6p . Inserting (IC.8I) and (IC.9I) in fjC.lj) and (IC.2[) . we obtain a 



system of two algebraic equations for the unknowns nE^ and J„ whose solution is 

^(o) _ aeajEp 

aj<y e + F 2 ' 



(CIO) 



A (0) = nE® - ijf = n a '^-f , (C.12) 
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The approximate electron current density (1C.11I) provides an approximate electron drift 
velocity vs. field, Vd(F) = Jn^ /n, whose maximum value is reached at 



Eg 

T 




Wo) 




v/ Ot e Otj , 



(C.13) 



*j 2/ (/?o) v 

in which we have used f l95l) to relate E to the lattice temperature l//3 = 2k B T /A for 
a Boltzmann local equilibrium. For a e — ctj — 1, (1C.10I) - flC.111) become the well-known 
values ( 142]) corresponding to the simple KSS-Poisson problem (J2J) - ((8]) with Boltzmann local 
equilibrium 18] provided r e = ^aj/a e . It is interesting to note that we have derived (IC.lOp 
and (IC.llj) for an unspecified general local equilibrium f 1Da , not just for the Boltzmann 
distribution. This means that this expression for the electron drift velocity is also valid 
at low temperatures, when the Fermi-Dirac distribution ( 189|) is a better description, and it 
justifies a posteriori the use of ( 1C.11I) to fit experimental results [301 ] . 



Remark CI. To leading order, E and J n in the right hand sides of ( IB. 21) and ( IB. 41) can be 
eliminated by using (IC.10I) and (1C.11I) . thereby obtaining 



jj* e p cos k cogh (uK) cos KdK 
J Q W e^ cosK cosh(uK) dK 



a e E 



aj + F 2 
aja e + F 2 



u 



e 13 smh(wr) 



« e (l - aj)E F 

Ct^Ole, + F 2 



f3 $ ^ eP cosK cosk{uK)dK ^ e 

Solving these two equations yield the functions (3(F) and u(F). In the case ctj 
yields u = and (IC.14j) becomes 

hifi) a e (l + F 2 )E 



W) 



a e + F 2 



(C.14) 
(C.15) 

i, (laT5|) 



(C.16) 



Equations (JHZD - (EHD yield 

A°) _ „ _ AD<*(o) Am) 

JO — II — Jo , Jo 

Ref[ m) = n£ (m) , Im/ 1 (m) 
The equations for /W and are 



0, for m = 1, 2, . . ., 



. 7"( m ) 



L/ (l) _ ylDa(l) = _ ^ 



+ sin k — — I 



1Dq(2) 
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+ sin k 



dx 

dx J dt 



(C.17) 
(C.18) 

(C.19) 
(C.20) 



and so on. The subscript m = 0,1 in the right hand side of these equations means that 
dF/dt and dn/dt are replaced by (J5 m o — J^i) and —dj^/dx, respectively. In these 
equations, the operator is defined by 



dn 

hu{k) = F—(k) + u(k). 
ok 



(C.21) 



The compatibility conditions flC.17p - flC.18p imply the following solvability conditions for the 
hierarchy flOT9l and flO20l : 

(L/^),- = 0, j =0,1. (C.22) 

Using the solvability conditions flC.22p for the linear hierarchy of equations, we can show 
that the reduced balance equations for n and F are obtained by inserting (1C.3I) in J n = — Im 
fi- 

oo 

J n = - y £ $ m Im/! (m) , J {m) = -Imf[ m) . (C.23) 

m=0 

We have already calculated = jf to be given by Eq. (IC.lip . To get a diffusive 

correction to this electron current density, we need to calculate Im/i . From (1C.19I) . (1C8H . 
fl09l and ffCT8]) . we obtain 



Refl 



(i) OjReri + Flmri 



a e ctj + F 2 
(i) a e lmri — F Reri 



a e ctj + F 2 



in which 



Thus we need to find 



0/(0) 



dt 



— sin k 



d_ n - ff Ofl 
dx 2i 



dx 



n = - 



r(0) 



in order to calculate (IC24I) . i.e., 



a. 



(i) 



Im 



(0) 



J 



Equation flCU2]) yields 



& 2 



— F 



Re 



(C.24) 
(C.25) 

(C.26) 

(C.27) 



(<» 



<:>/ 



IJ_T m f(°) 



0/i 



(0) 



dt 



a e aj + F 2 
+ n(J-jf) 



(aj-iF) 



dJ, 



(o) 



dx 

2acjF + i(a e aj - F 2 ) 
a e aj + F 2 



(C.28) 



(C.29) 
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The calculation of involves that of f2° a ^ . Using cos2/c = 1 — 2sin 2 /c, sin2fc = 
2 sin k cos k, integrating by parts and using (1B.2j) . flB.4|) from [Bj and ( 1C.14|) and (1C.15j) . 

we get 



n - Re ff _ a e nE 1 - (1 - a,)(l - uF) + F 2 







1 B nu 



aja e + F 2 

p2\~ 



a e nE (1 + F 2 )u - (1 - ay) [it + (1 + f3)F] 



For 1 — a,- 



(3 (3 
u = 0, we get Im/^ = and 
n-Re/ 2 B 



aj<y e + F 2 



(C.30) 
(C.31) 



a e nE 1 + F 2 



a e + F 2 ' 



(C.32) 



In this case, we obtain 



n - Refe 



(o) 



n 



— 2 lmf 2 



(o) 



1 + 4F 2 
nF 
1 + 4F 2 



a e E {l + F 2 
(a e + F 2 ) 
2a e F (l + F 2 
/3 (a e + F 2 ) 



2F 2 + 



1 - 



(C.33) 
(C.34) 



where $ is a function of F found by solving the equation (10.140 . 

Recapitulating, we have obtained the drift-diffusion equation flC4j) (Ampere's law) for 
F in which = is given by (ICTTTj) and is given by f l028|) - (l029|) and, in the 
particular case of a Boltzmann local equilibrium with ay = 1, by ( 1C.32I) - ( 1C.34j) . We have 

1 



+ Sim 



OF 






+ 


~dt 


df[ 0) 




dt 


0- 



a,a fi + F 2 



+ 5F 



a P 



E nF --A( n _ Re / 2 (0) ; 
2&r 1/2 ' 



2 9x 1/2 at 



j(t), 



(C.35) 



where n = 1 + dF/dx according to the Poisson equation (186]) . Note that the drift-diffusion 
equation ( 1C.35|) coincides with the drift-diffusion equation f)129p when we substitute hs = 
dfi/dt\ given by flC29[) (in which = O(S) has been neglected) and g s = ffl in f 1 1 2 9 j) 
with a e j = Sj ej j according to ( 1901) . Eq. ( 1C.35I) for almost elastic collisions becomes 



dF_ S 
~dt + F 2 + 5 2 7i7e 



j-t r\ r\ /• (0) 

7eE ° nF + 2 ^ Im /2 °* + 57eIm m 



0,5 



— — 7T— (n — Re f, q) — F Re — — 

2 dx y J2 ' SJ dt 



0,5 



(A 



(C.36) 
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where (IC.29I) should be inserted and is given by (IC.lj) and f )C2[) . In Equation ( 159"]) for 
the BO amplitude A, g = / 2 (0) + 5/ 2 (1) with f x = / 1>s + Ae~* e and / 1>s is given by (jCIQj) 
and (IC.llj) . Therefore in the case of almost elastic collisions, if the amplitude of the Bloch 
oscillations decays to zero, we are left with the above written drift-diffusion problem. 

Eq can be rewritten as 



Remark C2. Note that ([38]) with E 1D 

S^nEoiSjj - iF) 



fi,i 



F 2 + 5 2 7 ,7e 



+ 



+ 



F 2 + 5 2 ljle 
5 

F 2 + 5 2 ljle 



F + idj e d 
2 dx 



[n 



Re/ 2 , ) + 



{iF - 5 7i )Re^ — (F + z5 7e )Im 



6jj — iF d 

2 



Im / 2 ,c 



(C.37) 



5t v ,ey dt 

Our result for /i^ = df\ } s/dt means that we have approximated /^g by the first term in 
( 1C.37I) . The second and third terms in ( 1C.37I) correspond to ( 1C.27I ) which enter the 0(6) 
corrections (IC24[) and ( IC.25j) to the distribution function. Thus setting h s = dfi^s/dt 
corresponds to hs = (df[ ^ /dt)\ , with given by ( ]C.12j) . 
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